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This is a review of numerical applications of classical gluodynamics to heavy ion 
collisions. We recall some results from calculations of gluon production, discuss 
their implications for heavy ion phenomenology, and outline a strategy to calculate 
the number of quark pairs produced by these classical fields. 



1. Introduction 

To understand what is happening in relativistic heavy ion colliders, cur- 
rently RHIC in Brookhaven and in the future the LHC at CERN, one needs 
to understand not only hard probes, but also bulk particle production and 
thermalisation. It can be argued that the large phase space densities of 
partons in the small- a; wavefunction of the nuclei generate a hard enough 
momentum scale to allow weak coupling methods to be used. We shall 
discuss first some general ideas concerning weak coupling, classical field 
methods in this context. Then we will go on to discuss first gluon and then 
quark pair production in the McLerran-Venugopalan model. 



2. Relativistic heavy ion collisions 

We shall be interested in studying the case where two nuclei move at the 
speed of light along the x ± = 0-axes, i.e. at r = a . These nuclei then 
collide and leave behind them, at finite values of 77 and r, some matter 
which is then observed in detectors located in some region, varying between 
different experiments, around r\ — 0. We can divide the collision process in 
different stages: 



a The light cone coordinates are defined as = (t rt z)/\/2 and the proper time and 
spacetime rapidity as: r = \/t' 2 — z' 2 and n = t; In 3L_ . 



1 



February 2, 2008 0:50 Proceedings Trim Size: 9in x 6in 



lappi 



2 

(1) The initial condition at r = depends on the properties of the 
nuclear wavefunction at small x. 

(2) The thermal and chemical equilibration of the matter formed at 
r < t requires understanding of time dependent, nonequilibrium 
Quantum Field Theory. 

(3) The Quark Gluon Plasma, surviving for some fermis around To < 
t < lOfm. If the system reaches local thermal equilibrium, finite 
temperature field theory and relativistic hydrodynamics can be used 
to describe its behaviour. 

(4) Finally, for r > lOfm the system hadronises and decouples. 

The question of the thermalisation timescale To remains poorly under- 
stood. Hydrodynamical calculations have bees very successful in explaining 
the experimental observations, but their success depends on the assump- 
tion of a very early thermalisation time 1 . Most perturbative estimates, i.e. 
the bottom-up scenario 2 , generically produce quite a large thermalisation 
time, To > 3fm. On the other hand, one could argue that if the behaviour 
of the system is characterised by some quite large momentum scale, i.e. the 
saturation scale Q s <~ 1 . . . 2GeV, thermalisation should occur already at 
times t <~ 1/Q S ~ 0.2fm. It has been pointed out recently (see e.g. Ref. 3) 
that plasma instabilities could provide the rapid thermalisation that hy- 
drodynamical models require. In this context classical field models of the 
nuclear wavefunction and particle production can provide some insight into 
understanding the collision process. 

3. Saturation and the classical field model 

The general idea of parton saturation is that the the small- a; components of 
the nuclear wavefunction are dominated by a transverse momentum scale, 
the saturation scale Q s . The scale Q s is supposed to grow for decreasing x as 
Qs(x) <~ x~ x with A w 0.3 giving a good fit to HERA data on deep inelastic 
scattering 4 . Thus for small enough x or, equivalently, large enough y/s we 
have Q s 3> Aqcd and we can use weak coupling methods. The presence of 
the scale Q s is supposed to be caused by the color fields in the nuclear wave 
function becoming so strong that the gluonic interactions are dominated by 
the nonlinearities in the Yang-Mills Lagrangian. 

This idea can also be thought of as parton percolation. Let us attach, 
arguing by the Heisenberg uncertainty principle, to an individual parton 
with transverse momentum a transverse area ~ 1/Pt- One expects 
a qualitative change in the behaviour of the system at momentum scales 
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where the partons overlap and percolate the transverse plane: N/p\ ~ 
ttR\, where N is the number of partons and kR\ the nuclear transverse 
area. 

Because the color fields are strong, the occupation numbers of quantum 
states of the system are large, and it is natural ot use a classical field 
approximation. The saturation model has been dubbed the "Color Glass 
Condensate" and even called "a new form of matter". In theory the term 
"color glass condensate" refers to the small- a; wavefunction of the nucleus, 
characterised by the saturation scale Q s . In practice, most applications 
of these ideas to heavy ion collisions have so far not been classical field 
computations, but perturbative calculations with some phenomenological 
gluon fcx-distribution depending on Q 2 (x). 

3.1. Heavy ion collisions in the classical field model 

To study particle production in the classical field model we have to solve 
the Yang-Mills equations of motion: 

[D^F^] = r, (i) 

with a current of two infinitely Lorentz-contracted nuclei moving along the 
two light cones: 

J" = 5»+p {1) (x T )S(x-) + 5»-p {2) (x T )5(x+). (2) 

The model part of this approach comes when one must take some form for 
the charge density p{xt). The suggestion of McLerran and Venugopalan 5 
was to take a random stochastic color source with a Gaussian distribution: 

(p a (x T )p b (y T )) = g 2 p 2 S ab S 2 (x T y T ) (3) 

and then average all quantities calculated from p(xr) with this distribu- 
tion. An important development has been to consider a more general 
x-dependent probability distribution W x [p(a;T)] and derive a renormalisa- 
tion group equation for this distribution as a function of x, the JIMWLK 
equation 6 . 

3.2. Note on different saturation scales 

Different ways of understanding saturation have led to different definitions 
of the saturation scale in the literature, which can be a source of consider- 
able confusion. 
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The p sa t in the EKRT 7 model is a final state saturation scale, not a 
property of the nuclear wave function (initial state), and thus con- 
ceptually a bit different from the saturation scale we have discussed 
so far. 

We are using the old notation of of Krasnitz et. al. In their newer 
work 8 they define A s = g 2 ^. 

By calculating the gluon distribution in the McLerran-Venugopalan 
model one can relate the strangth of the color source to the satu- 
ration scale Q s defined by A. Mueller and Yu. Kovchegov (see e.g. 
Rcfs. 8, 9) b : 

«-^-(&) 

The saturation scale, or radius, Q s = l/i? s in the work of Golec- 
Biernat and Wusthoff 4 or Rummukainen and Weigert 10 is the same 
except with Ca replaced by Cf- c 

For a comment on the relation to the saturation scale used by E. 
Iancu et. al., see Ref. 11. 



4. Gluon production in the classical field model 



(1) 

An, = pure gauge 1 



(2) z_ 
Aa = pure gauge 2 



The solution of the Yang- 
Mills equations in regions (1) 
and (2) is an analytically 
known pure gauge field 12 , 
and gives the initial condition 
for the numerical solution in 
the forward light cone (3). 

The numerical method for solving the Yang-Mills equations in the forward 
light cone was developed by Krasnitz, Nara and Venugopalan 13,14 . The 




b Note the dependence on an infrared cutoff Aqce> that gives a large numerical uncer- 
tainty. 

c The techical reason for this is that they consider correlators of Wilson lines in the fun- 
damental representation, whereas the the gluon distribution that Mueller and Kovchegov 
consider involves the same correlator in the adjoint representation. 
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correct result for the transverse enegy was found in Ref. 15, (see also the 
erratum to Ref. 14 in Ref. 16). 

The numerical computation is done in the Hamiltonian formalism. Due 
to the boost invariance of the initial conditions the Yang-Mills equations 
can be dimcnsionally reduced to a 2+1 dimensional gauge theory with an 
adjoint scalar field. With the assumption of boost invariance one is explic- 
itly neglecting the longitudinal momenta of the gluons; a restriction that 
should be relaxed in future computations. In the Hamiltonian formalism 
one obtains directly the (transverse) energy. By decomposing the fields in 
Fourier modes one can also define a gluon multiplicity corresponding to the 
classical gauge fields. 

4.1. Numerical results 

Let us define the dimensionless ratios describing the energy and multiplicity: 

dE/d V dN/d V 

These dimensionless ratios depend, apart from the lattice spacing a, only 
on one dimensionless parameter characterising the field strength, g 2 fiR,A- 
As can be seen from Fig. 1, for strong enough fields (g 2 ^iRA > 50) both f E 
and /at are approximately independent of g 2 (J,RA and the lattice spacing. 
Figure 2 shows the energy as a function of time in different field components 
and the spectrum of the produced gluons. 
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Figure 1. Left: The dependence of /e and /jv on the field stregth paramater g 2 fiRA. 
Right: The dependence of f& and /jv on lattice spacing for different values of fi. 
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Figure 2. Left: the energy per unit rapidity is evenly distributed between different field 
components and almost constant after a very short time of the order of . This means 
that the 3+1 dimensional energy density decreases as e ~ r Right: The differential 
multiplicity for two different transverse lattice sizes. The curves labeled "KNV" are a 
fit to the numerical result of Ref. 14, and KNV "scaled" after these results have been 
corrected (see Ref. 16). 



4.2. Phenomenology 

At the level of the present discussion g 2 /i is still a free parameter that needs 
to be fixed from the experimental data. One can distinguish between two 
broad scenarios for relating the measured results to the calculated initial 
state quantities. 

Hydro scenario If the system thcrmalises fast, one can use ideal hydro- 
dynamics to follow its subsequent evolution. In this scenario entropy and 
thus multiplicity are approximately conserved, but the transverse energy 
decreases by a factor of <~ 3 due to pdF-work going down the beampipc. 
The corresponding value for the saturation scale is g 2 fi ~ 1.9GeV. 

Free streaming scenario Assuming that the gluons interact very 
weakly and do not thermalize, one can argue that longitudinal pressure 
is negligible and thus the transverse energy is conserved. Assuming energy 
conservation one gets a lower value for the saturation scale: </ 2 /z ~ 1.4GeV. 
This lower value means that the multiplicity must increase during thermal- 
isation or hadronisation of the system approximately by a factor of 2. 

At least in the first case a prediction for the multiplicity at central 
rapidities at the LHC is straightforward; one expects to see ( T^olS^r ) x 
1000 i=a 3000 particles per unit rapidity d . 



d Total, including neutral particles. At this level of approximation, N ch Ri ^N tot . 
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5. Quark pair production 

Given the classical fields corresponding to gluon production a natural ques- 
tion to ask is: how are quark- antiquark pairs produced by these color fields? 
Formally quark production is suppressed by a s and group theory factors 
compared to gluons, so in a first approximation we should be able to treat 
the quarks as a small perturbation and neglect their backreaction on the 
color fields. Heavy quark production is calculable already in perturbation 
theory 6 , but one can ask how much the strong color fields change the result. 
Understanding light quark production would address the question of chem- 
ical equilibration; turning the color glass condensate into a quark gluon 
plasma. 

The calculation of quark pair production, outlined in more detail in 
Rcf. 18, proceeds by solving the Dirac equation in the background color 
field of the two nuclei. This can be done analytically for QED . The initial 
condition for t — > — oo is a negative energy plane wave. Similarly to the 
QED case, one can find analytically the solution for the regions > 
0, :r+ < 0. These then give the initial condition at t = for a numerical 
solution of the Dirac equation in the forward light cone r > 0. To find the 
number of quark pairs one then projects the numerically calculated wave 
function to a positive energy plane wave at some sufficiently large time r. 

A major technical challenge in this calculation is the coordinate system. 
In order to include the hard sources of the color fields, the colliding nuclei, 
only in the initial condition of the numerical calculation, one wants to 
use the proper time r instead of the Minkowski time t. Unlike the gluon 
production case, where one was able to assume strict boost invariance, one 
now has a nontrivial correlation between the rapidities of the quark and the 
antiquark. Thus, although the background gauge field is boost invariant, 
one must solve the Dirac equation 3+1 dimensions. 

5.1. 1+1-d toy model 

The longitudinal direction is the one posing the most technical problems. 
To understand how to handle the longitudinal dimension we can construct 
a 1+1-dimensional toy model without the transverse dimensions. In 1+1- 



e In the weak field limit quark pair production from this classical field model reduces to 
a known result in fcx-factorised perturbation theory 17 . 

f The QED calculation, of interest for ultraperipheral collisions, is done e.g. by Baltz 
and McLcrran 19 and others. Baltz, Gelis, McLerran and Peshier 20 discuss the theory in 
more detail. 
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dimensions Dirac matrices are 2-dimensional. In the temporal gauge A T = 
that we have been using throughout the calculation there is only one 
component, A v , in the external gauge field. The mass m 2 s in the 1+1- 
dimensional model corresponds to the transverse mass k\ + m 2 of the full 
theory. 

The initial condition at r = involves a longitudinal momentum scale 
(the longitudinal momentum of the incoming antiquark) and thus we must 
use a dimensionful longitudinal variable to be able to represent the initial 
condition. We choose to take as our coordinates r and z. The Dirac 
equation becomes 

d T i> = Vr2 + Z2+7VZ (- 7 °7 3 5^ - mWV) - H°7 3 3^. (6) 

T T 

Because of the way the coefficients depend explicitly on the coordinates the 
discretisation of this equation is potentially very unstable and we have to 
use an explicit discretisation method. 
Let us take the background field as s : 

A v = cQ s T.h{Q s T). (7) 

For weak fields c « 1 we can compute the amplitude for quark pair 
production using the first order in perturbation theory (diagram (a) in 
Fig. 3). The result is a peak at 

2k+k- = (p + q) 2 = 2m 2 cB (1 + cosh(Ay)) = Q 2 S . (8) 

As can be seen from Fig. 4, for weak fields our numerical computation re- 
produces this perturbative peak. For stronger fields the numerical solution 
sums over all the diagrams in Fig. 3 and the position of the peak is shifted. 
In the full 3+1-dimensional case this peak is washed out by integration over 
the relative transverse momentum of the pair. 

A rapid back-of-the envelope estimate suggests that with transverse 
lattices of the order of 256 2 points (lattices from 128 2 to 256 2 were used in 
the gluon production computation) one could just manage to have a large 
enough lattice in the z-direction with the computers available to us. Then 
the computation should in principle be repeated for all the different values 
of pt (of the antiquark) . As this would be prohibitively expensive in terms 
of CPU-time one will, however, probably have to do with an interpolation 
from a reasonable amount of points on the pt lattice. 



g This is the correct time dependence of the perturbative solution to the gauge field 
corn's 12 . 
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Amplitude for different values of O 
cQ = 0.05 m s 



Amplitude for different values of Q 

cQ = m s 




Figure 4. Absulute value of the quark pair production amplitude for different values of 
the oscillation scale Q B . Left: weak fields, the peaks are at the location explained by 
Eq. (8). Right: strong fields. 



6. Conclusions 

Production of gluons has been computed numerically in the classical field 
model for heavy ion collisions, and the results are in rough agreement with 
RHIC observations, although not very conclusive due to the uncertainty in 
the numerical value of the saturation scale. A numerical calculation of the 
number of quark pairs produced from these classical fields is under way It is 
hoped that this computation will tell us something about the how a purely 
gluonic system can transform into a plasma of both quarks and gluons. 
Understanding kinetic thermalisation in this context will require treatment 
of the longitudinal dimension and implementing initial conditions from the 
JIMWLK equation. 
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